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Abstract 

The zero temperature properties of interacting 2 dimensional lattice bosons 
are investigated. We present Monte Carlo data for soft-core bosons that 
demonstrate the existence of a phase in which crystalline long-range order 
and off-diagonal long-range order (superfluidity) coexist. We comment on 
the difference between hard and soft-core bosons and compare our data to 
mean-field results that predict a larger coexistence region. Furthermore, we 
determine the critical exponents for the various phase transitions. 
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The possibility of a phase in dense Bose systems in which diagonal and off-diagonal 
long-range order (LRO) coexist has been the subject of discussion over the past 25 years [|T|. 
Normally bosons at zero temperature are either superfluid (with off-diagonal LRO) or solid 
(with diagonal LRO). However, for a finite range of the interactions between the bosons 
a coexistence phase was predicted within a mean- field approximation Experiments 
have been performed on A He, but no positive identification of this coexistence phase (often 
called supersolid) has yet been made. There are, however, strong hints towards such a phase 
||. On the theoretical side the discussion was restricted to the mean-field level. We are 
not aware of any more rigorous studies that identified a supersolid phase. In this letter we 
report on Monte Carlo simulations of soft-core lattice bosons in 2 dimensions that clearly 
demonstrate the existence of the supersolid phase beyond the mean-field approximation. 

Our lattice boson approach applies especially to artificially fabricated Josephson junction 
arrays (JJA) in which the bosons are Cooper pairs that tunnel between superconducting 
islands. In these arrays hopping, Coulomb interactions and the chemical potential may 
be tuned independently and therefore they constitute an ideal system for investigating a 
supersolid phase. Moreover, the Cooper pairs in JJA's are soft-core bosons that have a 
larger coexistence region than hard-core bosons as was demonstrated on the mean-field level 
by Roddick and Stroud M. 

The specific model we investigate is defined by the Hamiltonian for a JJA in the presence 
of an offset charge q that corresponds to a chemical potential for Cooper pairs || 

ij <ij> 

The number of excess Cooper pairs and the phase of the superconducting order parameter 
on island i are denoted by qi and (pi. Number and phase are conjugate variables that satisfy 
the commutation relation [qi,<f)j] = iSij. The average density of bosons may be varied by 
applying a gate voltage Vo, which enters in the chemical potential qo = CoVb/2e through the 
capacitance C between the islands and the ground ]7) . The Coulomb interaction Uy between 
Cooper pairs is determined by the inverse of the matrix that describes the distribution of 
capacitances between the islands. We take Uj to be short range, i.e. on-site and nearest 
neighbor interactions, Uq and U\, only. A natural stability condition is that the number of 
nearest neighbors times U\ has to be smaller than Uq. For Cooper pairs Uq is related to the 
charging energy E c hyU = 8E C = 4e 2 /C . If the Josephson coupling energy Ej dominates 
over the Coulomb interaction U^ the array will be superconducting at low temperatures. If 
on the other hand the Coulomb interaction dominates a Mott-insulator will be formed ||[7[ . 
The model defined by eq.(0) may be mapped onto a Bose-Hubbard model if the Josephson 
coupling term is identified with the hopping-term M. 

To gain understanding of the zero temperature properties of the model described by the 
Hamiltonian (|TJ), we first discuss the mean-field phase diagram following Roddick and Stroud 
0. The phase diagram is shown in figure 1 a) and b) for Ui/Uq=0.125 and 0.2 respectively. 
It is periodic in go with period 1 and symmetric around qo = ~. We discern four different 
phases: the superconducting phase (I), two incompressible Mott-insulating phases (II and 
III) and a compressible supersolid phase (IV). Phases I and IV have a nonzero superfluid 
density p s . Phases III and IV have nontrivial charge order ('checkerboard', see the inset to 
figure 1) and therefore a nonzero (7r,7r)-component of the static structure factor S n . Thus, 



2 



in the supersolid phase LRO (SV ^ 0) and off-diagonal LRO (p s ^ 0) coexist. 

The phase diagram for hard-core bosons was investigated in Refs. In that limit a 

supersolid phase is possible only in the presence of next nearest neighbor interactions. The 
difference with soft-core bosons is the lack of multiple occupation. Indeed, the expectation 
value for 2 soft-core bosons to be at the same site is nonzero in phase IV in figure 1 ||. 
We conclude that the possibility for bosons to hop over or past each other enhances the 
supersolid phase. 

The points marked a, (3 and 7 in figure 1 have particle-hole symmetry. This means 
that the cost in electrostatic energy is the same for adding or removing a boson. Point 
a and the phase boundary between phases I and II were investigated in refs. ||[7|]. Point 
a and the line separating phases I and II have a different dynamical critical exponent z. 
This exponent determines the space-time asymmetry. The correlation length in the time 
direction diverges like £ T ~ £ z , if £ is the correlation length in the space directions. Due to 
particle-hole symmetry the superconductor-insulator transition at point a has a dynamical 
critical exponent z=l. The transition is in the 3D XY universality class. For q ^0 the 
transition has z=2 and mean-field exponents apply. The same holds for point (3 and the line 
separating phases III and IV J7|. Motivated by these observations we expect also point 7 to 
have z=l, whereas for the transition at go 7^ § from phase I to IV we expect z=2. Below 
we show that this is consistent with our Monte Carlo data. The points marked 5 in figure 
1 have a first order transition, as the density jumps from in phase II to | in phase III. 

Since fluctuations around the mean-field solution are likely to be important in 2 dimen- 
sions one might wonder if the supersolid phase survives in an exact treatment. To investigate 
this question we performed Monte Carlo simulations of the model described by the Hamil- 
tonian flU). We follow closely the method used by S0rensen et al. Thus, we map our 2 



dimensional quantum model onto a 3 dimensional classical model of divergence-free current 
loops (we use the Villain form jll] for the cosine in eq.(|I]), see Refs. [0,0 for a derivation). 



The relevant quantity is then the partition function 



Yl ex p 

{J^=0,±l,±2,..} 




(2) 



where the sum is over divergence-free discrete current configurations that satisfy V M = 
(/i=x,y,r) and 5 < ij > equals 1 for nearest neighbors and is zero otherwise. The time- 
components of the currents correspond directly to the particle numbers, JJ = qi. The 
coupling constant K = 8fEj/Uo = fEj/Ec, where / depends on the time-lattice spacing. 
Here / is smaller than, but of the order of unity P,pT|. 

Using the standard Metropolis algorithm we generate configurations of currents in a 
system of size LxLxL T with periodic boundary conditions. The generation of configurations 
may be done canonical as well as grand-canonical. Here we work in the grand-canonical 
ensemble at fixed go m order to make contact to the phase diagrams in figure 1. As we 
are interested in a possible supersolid phase, the relevant quantities to measure are the 
superfluid density for off-diagonal LRO and the structure factor for diagonal LRO. 

In terms of the currents J M the superfluid density is 

P. = ^|lXf, (3) 
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It satisfies the finite size scaling relation |14| p s = L 2 d z p(bL 1 / u 5, L T /L Z ) with p a universal 



scaling function, b a nonuniversal scale factor, v the coherence-length critical exponent and 
5 = {K — K*)/K* the distance to the transition. At the critical point K = K*, 5 = and 
L z p s is a function of L T /L Z only. Thus, plots of L z p s vs. K will intersect at the transition 
if L T /L Z is kept constant. Furthermore, the data for L z p s plotted as a function of L l / y 8 for 
different system sizes should collapse onto one single curve. This allows the exponent v to 
be determined. 

A similar scaling relation holds for the structure factor [|15|] 



^ = -^£(-iy +i ^, (4) 

i.e. SV = L~ 2I3 I U S(U L l l u 5,L T / L z ), with the order parameter exponent (3. 

In the simulations we took Ui/Uq=0.2 in order to have a large coexistence phase. We 
performed simulations for constant qo= 0.5, 0.4 and and varied the coupling K. In the 
phase diagrams in figure 1 this corresponds to moving on horizontal lines through the phase 
transition(s). For q$= 0.5 and we simulated LxLx L systems, where L= 4, 6, 8, 10, 12, as 
suggested by particle-hole symmetry and z=l. Typically 100,000 sweeps through the lattice 
were needed for equilibration and the same amount for measurement. For qo= 0.4 we have 
z=2. In order to keep the ratio L T /L Z constant, we simulated L x L x L 2 /4 systems, where 
L= 6, 8, 10. For the largest system with L T = 25 we made up to 400,000 sweeps through 
the lattice for equilibration and the double for measurement. The results are summarized 
in figures 2-5 and table I. 

First we discuss our data for go=0.5. Figure 2 shows that there are two separate tran- 
sitions for diagonal and off-diagonal LRO with a coexistence region in between where both 
the superfluid density and the structure factor scale to a finite value in the thermodynamic 
limit. This demonstrates the coexistence of diagonal LRO and off-diagonal LRO for soft-core 
bosons with nearest neighbor interaction in 2 dimensions. Fluctuations reduced considerably 
the thickness of the supersolid phase compared to the mean-field result in figure 1 b. In fig- 
ures 3 a) and b) we plot again the same data around the critical points. In the neighborhood 
of the critical points the data fall onto a single curve when plotted as a function of L l l u 5. 
Table I shows that the exponent v is different for the two transitions. For the transition 
related to superfluidity (point (3 in figure 1) we find a value for v that is consistent with the 
3D XY universality class which has v «0.67. For the transition related to crystalline order 
(point 7) we find v «0.55 and (3 ~0.21. 

Also at q =0A we find two separate transitions that are the boundaries for the supersolid 
phase in between, see figure 4 and table I. As compared to go=0-5 both transitions are 
shifted to smaller values of the coupling constant K. This is consistent with the mean-field 
phase diagram. Again the two transitions have different critical exponents. The transition 
related to superfluidity (the line separating phases III and IV in figure 1) has v ^0.5 which is 
consistent with a mean-field transition in d+z=4 effective dimensions. The transition related 
to crystalline order (between phases I and IV) has an order-parameter exponent (3 RiO.25. 
This rules out a mean- field transition for diagonal LRO, although the transition is effectively 
4-dimensional. In the neighborhood of this transition, fluctuations of the x,y-components of 
the currents J induce long-range interactions in the time direction for the r-components of 
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the currents J @ . It is likely that these long-range interactions are a relevant perturbation 
and suppress the exponent f3. 

Finally the data for go=0 are shown in figure 5. Here there is only one phase transition, 
as the Mott-insulating lobes (phase II in figure 1) do not have any non-trivial crystalline 
order. Our data are consistent with a transition in the 3D XY universality class. 

In conclusion we have performed Monte Carlo simulations on soft-core lattice bosons 
in two dimensions that establish the existence of a supersolid phase in which diagonal and 
off-diagonal long-range order coexist. We estimated critical exponents as listed in table I. 
The mean-field phase diagram of ref. is qualitatively confirmed. However, our simulations 
indicate that the supersolid phase is smaller than one would deduce from mean-field theory. 
We suggest that the coexistence phase may be observed in two-dimensional systems such 
as Josephson junction arrays or thin 4 if e-films on suitable substrates. In these systems the 
possibility to vary the coupling constants as well as the chemical potential should make it 
possible to tune through the supersolid phase and see two sequential phase transitions. 

We acknowledge the help and suggestions of R. Fazio, C. Bruder, G. Schon and G.T. 
Zimanyi. This work is part of "Sonderforschungsbereich 195" which is supported by the 
Deutsche Forschungsgemeinschaft. 
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TABLES 





the transition for 


% 


off-diagonal LRO 


diagonal LRO 


0.5 


K* = 0.775 ± 0.005 


K* = 0.837 ±0.005 


(z = l) 


v = 0.65 ±0.08 


v = 0.55 ± 0.05, = 0.21 ± 0.04 


0.4 


K* = 0.645 ± 0.008 


K* = 0.749 ± 0.006 


(z = 2) 


v = 0.44 ±0.08 


v = 0.5 ± 0.1, = 0.25 ±0.10 


0.0 


K* = 0.843 ± 0.005 




(* = 1) 


^ = 0.61 ±0.08 





TABLE I. Critical couplings and exponents for the different transitions 



Figure Captions : 

Fig. 1 : Phase diagrams for soft-core bosons with on-site and nearest neighbor interaction, 
a) L r i/L r o=0.125. b) Ui/Uq=0.2. There are 4 different phases, I: superconductor, II: Mott- 
insulating, III: Mott-insulating with checkerboard charge-order and IV: supersolid. The 
points marked a, f3 and 7 have particle-hole symmetry. The points marked 5 have a first 
order transition. The inset to a) shows the checkerboard charge order. 

Fig. 2 : Data for Lp s and S n with —=0.78 vs. K at go=0.5. The curves cross at 
fC*=0.775 and 0.837 respectively. The region in between is the supersolid phase. 

Fig. 3 : Data for p s and S w at go=0-5 in the neighborhood of the two critical points, scaled 
as to collapse onto a single curve. The drawn lines are a low order polynomial fit to the 
data, a) Lp s vs. 5L l l v with i/=0.65. b) L°- 78 S n vs. 5L l l v with z/=0.55. 

Fig. 4 : Data for p s and S w at go=0.4. The drawn lines are a low order polynomial fit to 
the data, a) L 2 p s vs. K. The curves cross at fT*=0.645. b) L%S V vs. K with ^=1.0. The 
curves cross at K* =0.749. 



Fig. 5 : Data for p s at go=0. The drawn lines are a low order polynomial fit to the data. 
Lp s vs. K. The curves cross at K*=0.843. 
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